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Abstract: Masses associated with the topography, bathymetry, and its isostatic 
compensation are a dominant source of gravity field variations, especially at 
shorter wavelengths. On global scales the topographic/isostatic effects are also 
significant , except for the lowest harmonics. In practice, though, global effects 
need not be taken into account as such effects are included in the coefficients 
of the geopotential reference fields. On local scales, the short-wavelength gravi- 
ty variations due to the topography may, in rugged terrain, be an order of mag- 
nitude larger than other effects. In such cases, explicit or implicit terrain reduc- 
tion procedures are mandatory in order to obtain good prediction results. Such 
effects may be computed by space-domain integration or by FFT methods. 
Numerical examples are given in the paper for areas of the Canadian Rockies. 

In principle, good knowledge of the topographic densities is required to pro- 
duce the smoothest residual field. Densities may be determined from sample mea- 
surements or by gravimetric means, but both are somewhat troublesome methods 
in practice. The use of a standard density, e.g. , 2.67 g/cm 3 , may often yield 
satisfactory results and may be put within a consistent theoretical framework. 

The independence of density assumptions is the key point of the classical 
Molodensky approach to the geodetic boundary value problem. The Molodensky 
solutions take into account that land gravity field observations are done on a 
non-level surface. Molodensky's problem may be solved by integral expansions or 
more effective FFT methods, but the solution should not be intermixed with the 
use of terrain reductions. The methods are actually complimentary and may both 
be required in order to obtain the smoothest possible signal, least prone to alias- 
ing and other effects coming from sparse data coverage, typical of rugged topo- 
graphy. 


Introduction 


The two aspects of the role of the topography, namely the direct attraction 
of masses of the terrain and the uneven surface on which terrestrial measure- 
ments are made, are from a theoretical point of view completely different 
problems. The first assumes a density model, while the second - the Molodensky 
theory - in principle is free of any density assumptions. Using terrain reductions 
a computational smoothing of the gravity field is attempted, making interpolation 
and prediction from scattered data points more precise. Molodensky’s theory 
makes the classical geodetic boundary value problem solutions ’’correct” on the 
uneven topographic surface; applying the Molodensky correction terms to gravity 
implies no smoothing at al.. The methods are therefore complimentary and should 
be used together whenever feasible. 

Terrain reductions 


The terrain reductions may be classified under global or local models. On 
global scales, topographic-isostatic reductions must be used according to some 
idealized isostatic models. The simple model coming closest to geophysical rea- 
lity is the Vening Meinesz model, which is a modified Airy model taking into ac- 
count the elasticity of the crust, permitting short-wavelength loads to remain un- 
compensated. Some areas of the earth are, however, notably deviating from the 
simple models, as for example trench areas and midoceanic ridges. For a review 
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of computation of global topographic -isostatic effects, which may be done effici- 
ently by fast spherical harmonic expansions, see Rummel et al. (1988). 

The global topographic-isostatic reductions have only limited use in geodetic 
gravity field modelling, but are very useful for identifying anomalously compensa- 
ted areas. The real impact of terrain reductions come on the local scale, where 
they strongly diminish aliasing from undersampling of the rapidly varying, height- 
correlated gravity anomalies in rugged topography. 

On a local scale, isostatic effects need not be considered. Instead, a residual 
terrain model reduction ( RTM) , where only the topographic deviations frorn a smooth 
mean topography are considered, may be used with advantage. T e smoot mean 
height surface h r may be obtained from the topographic heights h by a moving 
average filter over suitable caps, say 1° in size. In this case, RTM-reduced gra 
vity anomalies will resemble isostatic anomalies (Forsberg, 1984). Removing t e 
complete effect of topography relative to a constant mean elevation level tor a 
given area, as often done in practice, may be considered a special case of the 
RTM reduction. The total computational removal of all topography down to t e 
geoid, i.e. the complete Bouguer reduction, is not useful in geodetic gravity field 
modelling because of very large indirect effects on the geoid. 

The general form of a terrain effect on any gravimetric quantity expressable 
as a linear functional L(T) of the anomalous potential T is of the general form 

h 

Lp(T m ) = G L p (-^)dxdydz, £= v (x-x p ) +(y-y p ) +(z-h p ) (1) 

E h r 

where E is the infinite x-y plane (planar approximation). Integrals like (l) may 
in practice be evaluated by prism integration or by expansions in co ™° “tions 
• permitting use of FFT techniques. For details see Forsberg (1984, 1985) or Side- 
ris (1985). Formally T m is the potential generated by the selected terrain mass 
model. When computations are done consistently in a remove-restore technique, 
i.e. modelling the reduced potential 

( 2 ) 


T = T - T 


m 


by subtracting terrain effects from input data, and restoring terrain effects in 
predictions, then in principle density p need not be known. However, meanmgfu 
results are only obtained when p is close to the real word values. 

The question of density 

Estimation of good insitu densities is often quite difficult. Only the surface 
of topography is available for sampling, and measurements of bulk densities on 
rock samples tend to show high variability even within the same geologica or 
mation. For sedimentary rocks, questions of porosity, water saturation, and com- 
paction present special problems. Typical density values encountered in practice 
range from below 2.0 g/cm 3 in moraine hills up to 3.3 g/cm in (rare) gabbroic 
intrusive areas. However, the standard density 2.67 g/cm’ represents a surprising 
good value in many cases (granities, gneisses, old sediments), and its use have 
been justified by many empirical investigations (Dobrin, 1976). Where significant 
ly lower density values need to be used (e.g., young sediments), topographic re 
lief is usually also lower and good density values are therefore less critical. 

For local applications a good alternative to density measurements is the estima- 
tion of p through studies of correlation of gravity with topographic heights ( the 
Nettleton method")- This method may be put within a consistent framework in 
least-squares collocation, estimating one or more density parameters alongside the 
gravity field modelling itself. For details, see Stinkel (19 ). 
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Molodensky’s problem 


The first order Molodensky solution to the geodetic boundary value problem 
consists of a series expansion of the form 


T, (P) = I 

n=0 2 hy J 


8 n (Q) 


[(xq-x p ) *(y 0 -y p ) ] 


2^ d V y Q 


(3) 


with 


g Q ( Q) = Ag(Q) 


g x (Q) 


-(hq hp) 2n 


Ag(x,y)- Ag(x n ,y n ) 

j ApS- dxdy 

[(xq-x) / +(yq-yn / 


(4) 

(5) 


For a review, see Moritz (1980) or Sideris (1987). Depending on the terrain roughness, 
higher-order terms may be quite significant. Their computation requires repeate 
applications of the harmonic continuation integral (5), which may be formulated 
as a sequence of convolutions and evaluated efficiently by FFT methods. For 

examples see Sideris (1987). , 

When used on free-air anomalies, formula (5) is known to be closely related 
to the classical terrain correction c (Moritz, 1980). The relationship comes from 
noting that in rugged topography free-air anomalies show a correlation with 
terrain height of the form Ag = Ag R +2aGph where Ag R is the simple Bouguer 
anomaly. Unfortunately the use of c ( requiring density information) rather than 
g (independent of p) have been the source of much confusion, and the practice 
is not recommended. It is thus preferable to evaluate (5) with terrain-reduced 
gravity data Ag c . The corresponding " reduced” Molodensky terms g n will be 
smaller, and the convergence of (3) improved. If terrain reductions are not used 
then, on the one hand, gravity data must be given densely enough to sufficiently 
sample even the shortest topographic - induced wavelengths, and, on the other 
hand, higher-order Molodensky series terms should be considered. Dense gravity data 
are hardly ever available in practice. 


An example: Kananaskis area, Canadian Rocky Mountains 

The Kananaskis area west of Calgary is a mountainous area with topography 
ranging from 1400 m to 3400 m. A number of astronomic deflections of the verti- 
cal and GPS-derived geoid undulations are available along the main valley area, 
in addition to gravity data space every 5 to 10 km in the surrounding region, 
and a dense digital terrain model. Results for a FFT gravity field modelling 
example are given below; for more results on terrain corrections for gravity and 
gradiometry in the same area, see Tziavos et al. (1988). 

For the terrain reduction of the available data, a 100 m x 100 m and 1 m 
x 1 km DTM was used. A smooth height surface with resolution of 70 km was 
generated by averaging. The statistics of the data, removing topographic effects 
relative to the mean height surface, is shown below. Predictions were not at- 
tempted without terrain reductions, as individual gravity anomaly values could chan 
ge up to lOOmgal close to the prediction points, depending on whether the obser- 
vation happened to be made on top of a mountain or at the bottom of a valley, n 
other words, the observed Ag— field is seriously undersampled without some kind 
of terrain reduction. 
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Kananaskis: Effect of the topography on observed gravity field data 


Data type 

*6 

(mgal) 

r, 

(arcsec) 

n 

(arcsec) 

C( m, relative to 

northernmost 

GPS-point) 

No. of 
points 

473 

15 


10 

Statistics 

mean std.dev. 

mean 

std.dev. 

mean std.dev. 

mean std.dev. 

Observed values 

20.5 

57.2 

1.37 

3.37 

.95 

6.28 

1.15 .46 

Terrain effects | 

-14.9 

47.2 

.86 

2.97 

-.22 

5.74 

.27 .21 

Reduced data 

35.4 

29.8 

.51 

1.60 

1.16 

2.15 

.88 .25 


For the FFT prediction, gravity anomalies were gridded on a 1.5' x 2. 5' grid 
in a 2.5 Q x 3° region. A similar 1.5' x2. 5' height grid was obtained by gridding 
gravity station heights for representing the uneven surface to which observations 
refer, in order to be able to use FFT techniques for computing the first-order 
Molodensky corrections g?. The procedure of gridding station heights rather than 
averaging the detailed DTM is preferable, because the height distribution of gra- 
vity stations does not necessarily follow the averaged topography; gravity stations 
tend to be located in valleys rather than mountaintops. 

The results of predictions with and without a spherical harmonic reference 
field (OSU86F to degree 180), and with and without the Molodensky terms, 
are shown below. Considering the rough topography, results are very satisfactory. 
The influence of the Molodensky term seems to be completely masked by other 
error sources, illustrating the high degree of smoothness locally provided by the 
terrain reductions. 


Kananaskis: Prediction results (observed minus predicted) for various solution types 


Solutior 

180x180 
reference field 

i with 

Molodensky 

terms 

t 

( arcsec) 

n 

(arcsec) 

C relative 
(m) 

mean 

std.dev. 

mean std.dev. 

mean 

KMo 

No 

No 

.92 

.65 

.31 .57 

.09 

.24 

No 

Yes 

.91 

.69 

.33 .60 

.09 

.25 

Yes 

No 

1.92 

.66 

i 

b* 

'•U 

‘s* 

00 

.09 

.17 
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